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Abstract 

Nicodemi and Prisco recently proposed a model for X-chromosome inacti- 
vation in mammals, explaining this phenomenon in terms of a spontaneous 
symmetry- breaking mechanism [Phys. Rev. Lett. 99 (2007), 108104]. 
Here we provide a mean-field version of their model. 

1 Introduction 

The nucleus of female mammals embryo cells contain two X chromosomes; one 
of these has to be inactivated in order to equalize the dosage of X genes prod- 
uct with respect to males, and this inactivation is crucial to survival [1, 2, 3]. 
Moreover, the X chromosome is conjectured to play a key role in the arising of 
certain types of cancer [4, 5]. 

The mechanism by which one of the two X chromosomes is inactivated is 
poorly understood, despite extensive work on this problem [1, 2, 3]; it is known 
that co-localization of chromosomes at the time X-chromosome inactivation 
(XCI) is established plays a key role [6, 7, 8]. 

Very recently, Nicodemi and Prisco [9] proposed a model to explain XCI 
on the basis of a spontaneous symmetry-breaking mechanism, related to the 
classical Ising model of Statistical Mechanics, in a dynamical stochastic model 
of the single-agent type. 

Their idea is that blocking factor (BF) molecules diffuse in the cell and 
can bind randomly to one or the other of the X-chromosomes; however, by a 
collective effect, the affinity to a chromosome increases when this already binds 
to other BF molecules. The net effect is that once one of the two chromosomes, 
by the effect of a fluctuation, binds to a larger number of BF molecules, it will 
bind to more and more of these - at a rate higher than the other chromosome 
(note that the BF molecules could have a greater probability to escape when 
more of them bind to the chromosome; this effect should be dominated by the 



growth in affinity in order for the collective effect devised by Nicodemi and 
Prisco to set in). In the whole, as a result of diffusive behavior, fluctuation and 
the affinity-enhancing collective effect, a net current will be established leading 
to BF concentration on one of the two X-chromosomes. 

The model can also account for the relevance of co-localization at the es- 
tablishment of XCI; in facts, if the two X-chromosomes are too far away the 
diffusive behavior will prevail, and the BF molecules escaping from the less 
populated X-chromosome will have a very low probability to reach the other 
one. 

The purpose of this short note is to provide a deterministic (simplified) ver- 
sion of the Nicodemi-Prisco symmetry-breaking XCI model, considering average 
exchanges of blocking factors between each of the two X-chromosomes and the 
cell fluid environment. Thus our model will amount to a mean- field version of 
the Nicodemi-Prisco model; we will formulate it first in terms of transition prob- 
abilities between states for a single BF molecule, and then pass to a formulation 
in terms of a first order system of ODEs, i.e. a Dynamical System [10, 11, 12], 
following a procedure which is standard in Physics or in Quantitative Biology 
[13,14]. 

Here we will consider the fluid environment as a homogeneous reservoir, i.e. 
overlook the effects connected with the spatial localization of the X-chromosomes; 
thus the model will be able to predict the symmetry breaking leading to inacti- 
vation of one of the X-chromosomes, but (in its simplest form, presented here) 
not its dependence on co-localization of the chromosomes. An extension of the 
model aiming at including co-localization effects will be developed in forthcom- 
ing work. 

Similarly here we will only deal with average quantities and flows; fluctua- 
tions can be included within this kind of modeling by promoting the ODEs we 
obtain to stochastic differential equations [15, 16, 17], as briefly discussed below. 

2 The model 

We will consider a fixed total amount (n units) of BF molecules; these can be 
in three mutually exclusive states: they can bind to one (state 1) or the other 
(state 2) of the two X-chromosomes, or be diffusing in the cell environment 
(state 3). We denote the number of units binding to the first X-chromosome at 
time t as X(t), that of units binding to the second X-chromosome at time t as 
Y(t), and that of units diffusing in the cell fluid at time t as Z(t). We will refer 
to x(t) and y(t) as the occupation numbers for the two X-chromosomes (or more 
precisely, for the Xic of the chromosomes) . Obviously, 

X(t)+Y(t) + Z(t) = N . (1) 
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2.1 Transition probabilities for a single BF molecule 

Let us consider the probabilities P(i — ► j, St) for transitions between these 
states - in particular, transition from state i to state j - for each BF molecule 
(referred to as a "particle" in the following) in a time interval of length St. In 
the following, we will use a simplified notation, and write 

Pij ■= P(i^j,St) . (2) 

We assume that for St sufficiently small, the probability a particle initially 
in state 3 will end up being in state 1 or 2 (i.e. bound to one of the two 
X-chromosomes) is 

pai a(X)St + o(St) ; P32 ~ a(Y)St + o{5t) . (3) 

That is, it will be proportional to the length of the time interval (but indepen- 
dent of time t itself) via a function a(x) or a(y) representing the affinity with 
the concerned X-chromosomc and depending on the number of particles already 
binding to it. 

Similarly, the probability that a particle initially binding to one of the X- 
chromosomes will escape from it in a time interval St will be given by 

p 13 ~ (3(X, Z)St + o(5t) ; p 23 ~ f3(Y, Z)St + o(St) . (4) 

That is, the escape probability will be proportional to the length of the time 
interval via a function of the number of particles binding to the concerned X- 
chromosome and of the number of particles fluctuating in the fluid environment 
(it is natural to assume these should enter only through their density z = Z/N). 

As for direct transitions from one to the other of the two X-chromosomes, 
these will be impossible on account of their non-zero space separation (the par- 
ticle will have to diffuse through the fluid to do that), so we will set 

P\2 = o(St) ; P21 = o(St) . (5) 

Needless to say, the probabilities of permanence in a given state are obtained 
simply requiring the conservation of particles, i.e. that the sum of probabili- 
ties for all transitions from a given state to all possible state is one (Markov 
property); thus 

pu = 1 - P12 - pis = 1 - P(X, Z)St + o(St) , 

P22 = 1 - P21 - P23 = 1 - P{Y A Z)St +0{St) , (6) 

P33 = 1 - Psi ~ P32 = 1 - [a(X) + a(Y)]St + o{5t) . 

2.2 Evolution equations for occupation numbers 

Let now {X(t), Y(t), Z(t)} be the average occupation numbers for states {1, 2, 3} 
at time t; we look at the average occupation numbers at time t + St. Note Z(t) 
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can be expressed using (1), hence we need only two equations. These will be 
given by 

X(t + 6t)= X(t)-p 13 X(t)+p 31 Z(t) , 

Y(t + St)= Y(t)-p 23 Y(t)+p 32 Z(t) . [<) 

Using the expressions given above for Pij - and omitting o(5t) terms - the 
(7) yield 

X(t + St) = X(t) + [a(X(t)) Z(t) - f3(X(t),Z(t)) X(t)]St 

Y(t + St) = Y(t) + [a(Y(t)) Z(t) - P(Y(t),Z(t)) Y{t)]5t . W 

Dividing by St and letting St — ► 0, eq.(8) yields a system of ODEs describing 
the evolution of the average occupation numbers for the three states (in which 
we make explicit the expression of Z{t) implied by the conservation law (1)): 

dX/dt= a(X)(N — X — Y) - (3(X,N — X — Y) X , 

dY/dt = a(Y) (N-X-Y) - /3(Y, N — X — Y)Y . [ ' 

Passing to consider the densities x = X/N and y = Y/N, and introducing 
the functions a and b defined through 

a(w) :=a(Nw) , b(w 1 ,w 2 ) ■= (3(Nw 1 , Nw 2 ) , (10) 

we get in the end the equations (symmetric under the exchange of x and y) 

dx/dt = a(x) (1 — x — y) — b(x,l — x — y)x 1 . . 

dy/dt = a(y) (1 - x - y) - b(y, 1 - x - y) y . 



2.3 Minimal model 

So far we have worked in completely general terms within the three-states de- 
scription of the system. In order to have a specific model, we should choose 
concrete forms for the functions o^iy) and (3{W) - and hence of the a(w) and 
b(w) - describing the probability of capture by and escape from chromosomes in 
function of the associated occupation numbers X, Y or densities x,y; this will 
yield a specific expression for the evolution equations (9) and (11). 

In their work, Nicodemi and Prisco [9] argued that the essential feature 
of their model for XCI is the collective phenomenon enhancing affinity and 
hence the probability of capture by X-chromosomes with a substantial number 
of binding particles. 

Thus the function a(W) should go to a finite limit ao for W — > (the BF 
have a non-zero affinity with the chromosomes even when no BF molecules are 
binding to chromosomes, and more generally when the collective behavior has 
not set in), and grow with W, i.e. a'(W) > 0. 

As for (3(W, Z), we should similarly have a function with a finite limit (3q for 
W — > 0; as for its dependence on W and Z, we assume it is only through the 
ratio X/Z of BF molecules binding to the chromosome and BF molecules in the 
fluid environment. 
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In the following, we choose to deal with a "minimal" model, i.e. consider a 
linear dependence on w for both these quantities; we set 

a(W)=a + a 1 W , 0(W) = b o + b 1 (W/Z); (12) 

this implies of course 

a(w) = a + aiNw , b(w) = b + b\{w/z) . (13) 

The parameters (do, a\, bo, b\) are all positive (we could always set one of 
these parameters equal to one by rescaling the time unit.) We also stress that 
in view of the discussion by Nicodemi and Prisco [9] , we should expect the coop- 
erative effects enhancing affinity should be predominant over those depressing 
the escape rate; this means we should expect a\ 3> |6i|. 

With these choices, the evolution equations for densities (11) read 

dx/dt= (ffl + cuNx) (1 - x - y) - [b + (&i/(l - N(x + y)))Nx] x , 
dy/dt = (a + ai Ny) (1 - x - y) - [b + (&i/(l - N(x + y)))Ny] y . 

(14) 

We are specially interested in the large N regime. In order to study this, it 
is convenient to rescale time (note that eqs.(14) become singular for N — > oo) 
and pass to consider as independent variable 

r := Nt . (15) 

We also set 

e := 1/N (16) 
(hence t — et); the equations (11) read then 

dx/dr = a\{l — x — y)x + e [a (l — x — y) — box] + e bix 2 (x + y — e)^ 1 , 
dy/dr= ai(l - x - y)y + e [a (l - x - y) - b y] + ebiy 2 {x + y - e)" 1 ; 

(17) 

at first order in e, we have 

dx/dr— oi(l — x — y)x + e [ao(l — x — y) — (bo — b\x/{x + y))x] , 
dy/dr= cn(l - x - y)y + e[a (l-x-y) - (b - hy/(x + y))y] . 

(18) 



3 Different time-scales for the dynamics 

We want now to discuss some point concerning the qualitative behavior of the 
dynamical system defined by (14); we should consider it in the region (invariant 
under the dynamics) of R 2 delimited by the coordinate axes and by the line 
x + y = 1 . 

As noted above, the symmetric fixed point E3 is a saddle point for the 
dynamics. It is easy to see that its (invariant) stable manifold [10, 11, 12] is 
just the line y = x. 
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Figure 1: Evolution of x(t) (upper curves) and y(t) (lower curves) as a function 
of the rescaled time r, obtained numerically integrating the equations (17). 
The parameters defining the model have been set as follows: do = 1, o-i = 0.01, 
b a = 0.5, 6i = 0.1; while for N we have used N = 10000, giving s = 0.0001. 
Initial conditions are x(0) = 0.011, y(0) = 0.01. (a): the initial expansion (this 
plot shows the evolution for < r < 10 3 ); note the two densities grow together, 
(b): on a longer timescale one observes how the slow phase sets in and the 
densities remain nearly constant over adong time (this plot shows the evolution 
for r < 5 • 10 4 ). (c): on a still longer timescale, the transition from the saddle 
point E 3 to the stable equilibrium - in this case, E 2 - is clearly visible and 
happens at an intermediate speed (this plot shows the evolution for r < 2 • 10 6 ). 
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Figure 2: Evolution of £(t) — x{t) + y{t) (dotted curves) and -q(t) = x{t) — 
y{t) (solid curves) as a function of the rescaled time r, obtained numerically 
integrating the equations (17). Parameters and initial conditions as in Figure 
1. (a): in the initial expansion £ grow abruptly, while ij remains near to zero 
(this plot shows the evolution for < r < 10 3 ). (b): in the slow phase both £ 
and r\ remain nearly constant over a long time (this plot shows the evolution for 
r < 5 • 10 4 ). (c): in the intermediate speed phase, £ remains nearly constant, 
while rj undergoes a relatively fast gr<7wth (this plot shows the evolution for 
r < 2 • 10 6 ). 
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Figure 3: Evolution of £'(t) (dotted curves) and rf{r) (solid curves) resulting 
from the numerical integration of the equations (17). Parameters and initial 
conditions as in Figure 1. (a): in the initial expansion £'(r) reaches values of 
order 2.5 • 1CP 3 , while ?/(t) remains smaller than 7 • 10~ 5 . (b): in the slow 
phase both speeds are very small, but £' < 10~ 8 while rj ~ 10~ 7 . (c): in the 
intermediate speed phase, £' remains extremely small (in this case we actually 
have £' ~ 10~ 9 ), while ?/ slightly grows again; in this case v( ~ 10~ 6 . 



Let us consider the situation where at first very few BF molecules are binding 
to either one of the X-chromosomes, i.e. x(0) ~0k 2/(0); note that as x(0) « 
y(0), the initial data are close to the stable manifold for the saddle point E3. 

In a first stage, both x(t) and x(t) will grow exponentially in t, with exponent 
01; this is easily seen from (14) (assuming i< 1, j;< 1). More precisely, near 
the origin one has dx/dr w aix, dy/dr w aiy; it follows that in this first phase, 
x(i) i=a 2/(f), i.e. the dynamics remains in a neighborhood of the symmetric line 
x = y. 

This fast growth will go on until the system gets near to £3 (note this also 
means 1 — (x + y) = 0(e)); at this point, a slower evolution gets in. Once it 
reaches the vicinity of £3, the dynamics is dominated by the expansion rate 
near £3, i.e. by the positive eigenvalue A 2 (i?3) w (bi/2)e, and evolution takes 
place at a speed of order e. 

However, at some point (after a time of order 1 je) the system gets far enough 
from £3 to have again a dynamics non fully described by the linearized system 
identified by M3, and a somewhat faster - albeit not as fast as in the first phase 
- evolution can take place. 

This will lead the system near one of the non-symmetric equilibria (depend- 
ing on the initial data), where again the dynamics is dominated by the eigenvalue 
which is smaller in modulus; this is A2(£i,2) = i-e. the system gets again a 
speed of order e, and thus approaches the final equilibrium, as already discussed, 
with a rate of order e. 

Summarizing, the evolution of our system will have (assuming initial condi- 
tions are near to zero) four distinct phases: 

(1) A phase of rapid growth, exponential with expansion rate of order a\; 

(2) A slow phase (speed of order e) spent in the vicinity of the saddle point £3; 

(3) An intermediate (moderate) speed phase, in which the system travels from 
a neighborhood of £3 to a neighborhood of either ones of the asymmetric stable 
equilibria £ 12 with a speed substantially higher than e; 

(4) A new slow phase, in which the system approaches exponentially the stable 
equilibrium at a rate of order e. 

This behavior is clearly shown in Figure 1, where we plot the numerical 
solution to the full equations for x(t) and y(r) (that is, without truncation to 
first order in e) on different time scales. 

In Figure 2, we look at the evolution of the quantities £(i) := x(t) + y(t) and 
rj(t) := x(t) — y(t), representing the total fraction of BF molecules binding to 
either one of the X-chromosomes and the symmetry breaking measure respec- 
tively. This again shows clearly that the overall evolution of the system goes 
through the different phases enumerated above. 

Confirmation is also obtained by Figure 3, where we plot d^/dr and drj/dr 
on different timescales. 
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4 Fluctuations 



The model equations (17) or the limit ones (18) only deal with evolution of 
average quantities, i.e. do not take into account in any way fluctuations. Our 
present task is to go beyond this level of description, and take into account 
fluctuations around the average dynamics. 

As a first approximation, these can be taken into account by transforming 
the equations (17) or (18) into stochastic differential equations [15, 16, 17] by 
adding a noise term (see the Appendix for details). In other words, if the 
deterministic equations are dx/dr — f(x,y), dy/dr = g(x,y), we will write 



with lu(t) a standard Wiener process and 7 a parameter describing the size of 
fluctuations. Note that (as these are fluctuations for an average) 7 will depend 
on N and hence on e; more precisely, 7 ~ 70 \fe. This dependence - and the 
estimate 70 < 7o = ~ can be obtained by a standard procedure. 

At first order in e, we have 



When l — x—y = 0(e), the noise term could - depending on the relation between 
e, ai and 70 - dominate the dynamics. That is, in the slow phase near the saddle 
point (see the discussion in section 3) random motion due to fluctuations could 
be dominant: if 7 is large enough, the time spent in a neighborhood B of E3 
will be determined by the escape time needed to exit B due to fluctuations. 
This will produce an acceleration of the dynamics with respect to the purely 
deterministic one - and possibly could cause the system to cross the line x = y 
and end up near the other stable equilibrium (this happens in the numerical 
simulation plotted in figure 4). On the other hand, if 7 is small enough then 
the fluctuations will cause just a twiggling of trajectories around deterministic 
ones. 

Obviously dominance of fluctuation in the dynamics will take place once the 
system gets near - or very near if 7 is small - to the stable equilibrium point E\ 
or E2. There will be a metastable equilibrium distribution around E\^ and - 
over extremely long times - large fluctuations could drive the system near to the 
other stable equilibrium; the equilibrium distribution and the time needed for 
large fluctuations to appear can be estimated via standard stochastic processes 
techniques [17]. 

In Figure 4 we show a numerical simulation of (the stochastic process defined 
by) the stochastic differential equations (20). 



dx = f(x,y)dr + ^dus^) , 
dy = g(x,y)dr + jduj(r) , 



(19) 



dx 



dy 



[ai (1 — x — y)x + e[ao(l — x — y) 

-(b - hx/(x + y))x}} dr + v / £7o^(t) , 

[ai(l -x-y)y + e[a (l -x-y) 

-(bo - hy/(x + y))y]} dr + y/e^ du(T) . 



(20) 
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Figure 4: Simulation of the stochastic process defined by (20); the parameters 
of the model have been chosen as in the previous Figures (see caption to Figure 
1 for their values), and we have set 7 = 10~ 4 . In this case the fluctuation do not 
seriously modify the dynamics, and the discussion based on the deterministic 
model confirms its validity. 

5 Conclusions 

We have provided a mean-field formulation of the model recently proposed by 
Nicodemi and Prisco [9] for the X-chromosome inactivation in mammals [2, 1, 



Like their model, the version proposed here explains X-chromosome inac- 
tivation as the result of a spontaneous symmetry breaking in the dynamics of 
blocking factors molecules binding to the X-chromosome inactivation centers 
(Xic). While their approach was based on a single-agent model, in our version 
the considered equations deal directly with average quantities. 

After discussing a general formulation of the model, we considered a minimal 
version, in which the affinity grows linearly with the number of BF molecules al- 
ready binding to the Xic and the rate of escape of BF molecules increases linearly 
on the ratio between density of BF molecules binding to the X-chromosomes and 
of BF molecules floating in the fluid environment with the same number. Our 
model contains four positive parameters, related to affinity and escape proba- 
bility for the BF molecules. 

We showed that when b\ < bo, non-symmetric equilibria exist; under the 
same condition they arc stable (while symmetric equilibria are unstable). Thus 
in such case we have a spontaneous symmetry breaking, i.e. non-symmetric 
equilibria describe the asymptotic state of the system. We were also able to 
estimate the rate of approach to equilibrium as a function of the parameters 
describing the system and the total BF population. 

A qualitative discussion showed that when the initial situation corresponds 
to most of the BF molecules floating in the fluid, one should expect the dynamics 



3, 6, 7, 8]. 
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to undergo different phases. We ran some numerical simulations (one of these 
is shown in the Figures) confirming the conclusions reached by the qualitative 
discussion mentioned above. We also briefly discussed how fluctuations can be 
taken into account within our scheme. 

Finally, we would like to stress that our model differs from that of Nicodemi 
and Prisco only in the mathematical description it considers, while the physical 
basis is the same. Thus, we are just providing an alternative description - in 
terms of dynamical systems and ODEs rather than of stochastic processes - for 
the Physics described by Nicodemi and Prisco in their paper [9]. 

This mathematically simpler description allows to make more detailed pre- 
dictions, in particular about different timescales in the dynamics; these were 
fully confirmed by numerical simulations. 

Our model, like the one of Nicodemi and Prisco [9] , did not take into account 
co-localization; work is in progress to take this into account within the present 
description of X-inactivation. 

Appendix 

In this Appendix we will briefly discuss the passage from the physical model of 
section 2 to the stochastic differential equations (19) and (20). 

In order to take into account full detail of fluctuations in our model, we 
should consider the flows $y of particles from state i to state j. Denoting by 
Vij (5t) the number of particles passing from state i to state j in a time interval 
of length det, we have of course 

®ij = v ij - v ji • ( 21 ) 

The transition probabilities pij for a single particle are given in section 2.1; the 
Vij will thus be Poisson distributed with parameter Ay = PijUi, where n, is the 
population of state i at the beginning of the time interval. 

Needless to say, if we look just at expectation values we obtain again (8); 
however, our present task is to go beyond this level of description, and take into 
account first and higher order momenta of the Poisson distribution - that is, in 
physical terms, fluctuations around the average dynamics. 

Let us fix our attention, for the sake of concreteness, on the population X (or 
density x = X/N) of the state 1 and its variation SX (or 5x) in a time interval 
of length St; we also write Z and z rather than expressing these quantities in 
terms of (X, Y) and (x, y), for ease of notation. 

We obviously have deX = z/31 — z/13. Recall that 1/31 is Poisson distributed 
with parameter A = a(X)Z5t (i.e. V(y^\ = k) = \ k e~ x /k\), and 1/13 is also 
Poisson distributed but with parameter fi = (3(X,Z)X. Thus $31 is Poisson 
with mean A — fx and variance a\ x = (A + Thus for X we have the stochastic 
differential equation dX = (A — ii)dt + adu>(t); in the case of our minimal model 
this reads 

dX = [(a +a 1 X)Z-(bo+b 1 X/Z)X]dt+y/(a + a x X)Z - (b + b 1 X/Z)Xduj(t) . 

(22) 
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Passing to density variables this reads dx = N 1 [(ao + a\Nx)Nz — (bo + 
bix/z)Nx]dt + iV _1 [(ao + a x Nx) - (b + b 1 x/z)Nx] 1 / 2 du>(t); we should also 
pass to the rescaled time variable r = Nt, so that dt — (1/N)(dr) and dw(t) = 
(l/s/N)du{T). 

In this way we get, using again TV -1 = e, 

dx = [aixz+e(a z — (b +bix/z)x]dT+y / 'aixz + e(aoz — (bo + b\x / z)x\fe <1uj(t) . 

(23) 

We thus have a noise of intensity 7 = 70 \fe\ note that at order zero in e, 
70 = yja\xz\ as the sum of x and z cannot be higher than one, this can be 
estimated by 70 < 7o = \fd\j1. 
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